import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import yfinance as yf
from scipy import stats
import warnings
warnings.filterwarnings('ignore')

plt.style.use('seaborn-v0_8-whitegrid')
plt.rcParams['figure.figsize'] = (12, 6)
np.random.seed(42)

print('Libraries loaded successfully!')

# Download stock data
tickers = ['SPY', 'AAPL', 'MSFT', 'GLD', 'TLT']
start_date = '2019-01-01'
end_date = '2024-01-01'

print(f'Downloading data for {tickers}...')
data = yf.download(tickers, start=start_date, end=end_date, progress=False)

# Handle different yfinance column structures
if isinstance(data.columns, pd.MultiIndex):
    if 'Adj Close' in data.columns.get_level_values(0):
        prices = data['Adj Close']
    elif 'Close' in data.columns.get_level_values(0):
        prices = data['Close']
    else:
        prices = data.xs('Close', axis=1, level=1) if 'Close' in data.columns.get_level_values(1) else data.iloc[:, :len(tickers)]
else:
    prices = data['Adj Close'] if 'Adj Close' in data.columns else data['Close']

returns = prices.pct_change().dropna()

print(f'Data range: {prices.index.min().date()} to {prices.index.max().date()}')
print(f'Trading days: {len(prices)}')
prices.tail()

# Let's look at SPY returns
spy_returns = returns['SPY']

# Calculate mean and median
mean_ret = spy_returns.mean()
median_ret = spy_returns.median()

print('=== SPY Daily Returns ===')
print(f'Mean:   {mean_ret:.6f} (annualized: {mean_ret*252:.2%})')
print(f'Median: {median_ret:.6f}')
print(f'\nDifference: {(mean_ret - median_ret):.6f}')

# Calculate volatility for all assets
print('=== Annualized Volatility ===')
print('(Higher = more risky)\n')

for ticker in tickers:
    daily_std = returns[ticker].std()
    annual_vol = daily_std * np.sqrt(252)
    print(f'{ticker}: {annual_vol:.1%}')

# Calculate skewness and kurtosis
print('=== Return Distribution Shape ===')
print(f'{"Asset":<6} {"Skewness":>10} {"Kurtosis":>10}')
print('-' * 28)

for ticker in tickers:
    skew = returns[ticker].skew()
    kurt = returns[ticker].kurtosis()
    print(f'{ticker:<6} {skew:>10.2f} {kurt:>10.2f}')

# Visualize SPY return distribution
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Histogram with normal overlay
ax1 = axes[0]
ax1.hist(spy_returns, bins=50, density=True, alpha=0.7, 
         color='steelblue', edgecolor='white', label='Actual SPY')
x = np.linspace(spy_returns.min(), spy_returns.max(), 100)
ax1.plot(x, stats.norm.pdf(x, mean_ret, spy_returns.std()), 
         'r-', lw=2, label='Normal Distribution')
ax1.axvline(x=0, color='gray', linestyle='--', alpha=0.5)
ax1.set_xlabel('Daily Return')
ax1.set_ylabel('Density')
ax1.set_title('SPY Returns vs Normal Distribution')
ax1.legend()

# QQ-plot
ax2 = axes[1]
stats.probplot(spy_returns, dist='norm', plot=ax2)
ax2.set_title('Q-Q Plot: Are Returns Normal?')

plt.tight_layout()
plt.show()

# How often do extreme events occur under a Normal distribution?
print('=== Probability of Extreme Events (Normal Distribution) ===')
print()

spy_std = spy_returns.std()
spy_mean = spy_returns.mean()

for sigma in [2, 3, 4, 5]:
    prob = 2 * (1 - stats.norm.cdf(sigma))  # Two-tailed
    expected_days = int(1 / prob) if prob > 0 else float('inf')
    print(f'{sigma}-sigma event: {prob:.2e} probability')
    print(f'   Expected once every {expected_days:,} trading days ({expected_days/252:.0f} years)')
    print()

# Count actual extreme events in SPY
print('=== Actual vs Expected Extreme Events in SPY ===')
print()

for sigma in [2, 3, 4, 5]:
    threshold = sigma * spy_std
    extreme_days = spy_returns[abs(spy_returns - spy_mean) > threshold]
    count = len(extreme_days)
    
    normal_prob = 2 * (1 - stats.norm.cdf(sigma))
    expected_count = normal_prob * len(spy_returns)
    
    print(f'{sigma}-sigma events: Actual={count}, Expected={expected_count:.1f}, Ratio={count/max(expected_count, 0.1):.1f}x')

# Compare Normal vs t-distribution
fig, ax = plt.subplots(figsize=(12, 6))

x = np.linspace(-5, 5, 1000)

ax.plot(x, stats.norm.pdf(x), 'b-', lw=2, label='Normal')
ax.plot(x, stats.t.pdf(x, df=3), 'r-', lw=2, label='t (df=3)')
ax.plot(x, stats.t.pdf(x, df=5), 'g-', lw=2, label='t (df=5)')
ax.plot(x, stats.t.pdf(x, df=10), 'orange', lw=2, label='t (df=10)')

ax.set_xlabel('Standard Deviations from Mean')
ax.set_ylabel('Probability Density')
ax.set_title('Normal vs Student\'s t-Distribution')
ax.legend()
ax.set_xlim(-5, 5)

plt.tight_layout()
plt.show()

print('Notice: Lower degrees of freedom = fatter tails = more extreme events')

# Fit t-distribution to SPY returns
df_fit, loc_fit, scale_fit = stats.t.fit(spy_returns)

print('=== Fitted t-Distribution Parameters ===')
print(f'Degrees of freedom: {df_fit:.2f}')
print(f'Location (mean):    {loc_fit:.6f}')
print(f'Scale (std):        {scale_fit:.6f}')
print(f'\nInterpretation: df={df_fit:.1f} confirms fat tails in SPY returns')

# Test for normality
print('=== Normality Tests ===')
print()

for ticker in tickers:
    ret = returns[ticker]
    jb_stat, jb_pval = stats.jarque_bera(ret)
    is_normal = 'YES' if jb_pval > 0.05 else 'NO'
    print(f'{ticker}: Jarque-Bera p-value={jb_pval:.2e}, Normal? {is_normal}')

# Test if SPY has statistically significant positive returns
print('=== One-Sample t-Test: Is SPY Mean Return Zero? ===')
print()

t_stat, p_value = stats.ttest_1samp(spy_returns, 0)

print(f'Sample mean:     {spy_returns.mean():.6f}')
print(f'Sample size:     {len(spy_returns)}')
print(f't-statistic:     {t_stat:.4f}')
print(f'p-value:         {p_value:.4f}')
print()

if p_value < 0.05:
    print('Result: REJECT null hypothesis - SPY has significant non-zero returns')
else:
    print('Result: CANNOT reject null hypothesis')

# Compare AAPL vs SPY
print('=== Two-Sample t-Test: AAPL vs SPY ===')
print()

aapl_ret = returns['AAPL']
spy_ret = returns['SPY']

t_stat, p_value = stats.ttest_ind(aapl_ret, spy_ret)

print(f'AAPL mean: {aapl_ret.mean()*252:.2%} annualized')
print(f'SPY mean:  {spy_ret.mean()*252:.2%} annualized')
print(f't-statistic: {t_stat:.4f}')
print(f'p-value:     {p_value:.4f}')
print()

if p_value < 0.05:
    print('Result: Significant difference between AAPL and SPY')
else:
    print('Result: No significant difference')

# Calculate correlation matrix
corr_matrix = returns.corr()

print('=== Correlation Matrix ===')
print(corr_matrix.round(3))

# Visualize correlations
fig, ax = plt.subplots(figsize=(10, 8))

im = ax.imshow(corr_matrix, cmap='RdBu_r', vmin=-1, vmax=1)

ax.set_xticks(range(len(tickers)))
ax.set_yticks(range(len(tickers)))
ax.set_xticklabels(tickers)
ax.set_yticklabels(tickers)

for i in range(len(tickers)):
    for j in range(len(tickers)):
        ax.text(j, i, f'{corr_matrix.iloc[i, j]:.2f}',
                ha='center', va='center', color='black', fontsize=12)

ax.set_title('Asset Correlation Matrix', fontsize=14, fontweight='bold')
plt.colorbar(im, label='Correlation')
plt.tight_layout()
plt.show()

# Demonstrate diversification
print('=== Diversification Benefit ===')
print()

spy_vol = returns['SPY'].std() * np.sqrt(252)
tlt_vol = returns['TLT'].std() * np.sqrt(252)
correlation = returns['SPY'].corr(returns['TLT'])

print(f'SPY volatility: {spy_vol:.1%}')
print(f'TLT volatility: {tlt_vol:.1%}')
print(f'Correlation:    {correlation:.3f}')
print()

# 50/50 portfolio
weighted_avg_vol = 0.5 * spy_vol + 0.5 * tlt_vol
portfolio_returns = 0.5 * returns['SPY'] + 0.5 * returns['TLT']
actual_portfolio_vol = portfolio_returns.std() * np.sqrt(252)

print(f'50/50 Portfolio:')
print(f'  If correlation=1: {weighted_avg_vol:.1%}')
print(f'  Actual:           {actual_portfolio_vol:.1%}')
print(f'  Benefit:          {weighted_avg_vol - actual_portfolio_vol:.1%} reduction')

# Rolling correlations
window = 60

rolling_corr_spy_tlt = returns['SPY'].rolling(window).corr(returns['TLT'])
rolling_corr_spy_gld = returns['SPY'].rolling(window).corr(returns['GLD'])

fig, ax = plt.subplots(figsize=(14, 6))

ax.plot(rolling_corr_spy_tlt.index, rolling_corr_spy_tlt, label='SPY-TLT', linewidth=1.5)
ax.plot(rolling_corr_spy_gld.index, rolling_corr_spy_gld, label='SPY-GLD', linewidth=1.5, alpha=0.8)

ax.axhline(y=0, color='gray', linestyle='--', alpha=0.5)
ax.set_xlabel('Date')
ax.set_ylabel('60-Day Rolling Correlation')
ax.set_title('Rolling Correlations Over Time')
ax.legend()
ax.set_ylim(-1, 1)

plt.tight_layout()
plt.show()

# Module Project: Your implementation here

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import yfinance as yf
from scipy import stats
import warnings
warnings.filterwarnings('ignore')

plt.style.use('seaborn-v0_8-whitegrid')
plt.rcParams['figure.figsize'] = (12, 6)
np.random.seed(42)

print('Libraries loaded successfully!')

# Download stock data
tickers = ['SPY', 'AAPL', 'MSFT', 'GLD', 'TLT']
start_date = '2019-01-01'
end_date = '2024-01-01'

print(f'Downloading data for {tickers}...')
data = yf.download(tickers, start=start_date, end=end_date, progress=False)

# Handle different yfinance column structures
if isinstance(data.columns, pd.MultiIndex):
    if 'Adj Close' in data.columns.get_level_values(0):
        prices = data['Adj Close']
    elif 'Close' in data.columns.get_level_values(0):
        prices = data['Close']
    else:
        prices = data.xs('Close', axis=1, level=1) if 'Close' in data.columns.get_level_values(1) else data.iloc[:, :len(tickers)]
else:
    prices = data['Adj Close'] if 'Adj Close' in data.columns else data['Close']

# Calculate returns
simple_returns = prices.pct_change().dropna()
log_returns = np.log(prices / prices.shift(1)).dropna()

print(f'Data range: {prices.index.min().date()} to {prices.index.max().date()}')
print(f'Trading days: {len(prices)}')
prices.tail()

# Calculate simple returns
print('=== Simple (Arithmetic) Returns ===')
print()
print('First 5 days of SPY simple returns:')
print(simple_returns['SPY'].head())
print(f'\nSum of all simple returns: {simple_returns["SPY"].sum():.2%}')

# Compare sum of returns vs actual total return
spy_prices = prices['SPY']

# Actual total return
actual_total = (spy_prices.iloc[-1] / spy_prices.iloc[0]) - 1

# Sum of simple returns (WRONG approach)
sum_of_returns = simple_returns['SPY'].sum()

# Compounded returns (CORRECT approach)
compounded = (1 + simple_returns['SPY']).prod() - 1

print('=== Total Return Calculation ===')
print()
print(f'Actual total return:     {actual_total:.2%}')
print(f'Sum of simple returns:   {sum_of_returns:.2%}  ← WRONG!')
print(f'Compounded returns:      {compounded:.2%}  ← CORRECT!')
print()
print('Lesson: Never sum simple returns to get total return!')

# Log returns ARE additive across time!
actual_log_total = np.log(spy_prices.iloc[-1] / spy_prices.iloc[0])
sum_of_log = log_returns['SPY'].sum()

print('=== Log Returns: Time Additivity ===')
print()
print(f'Actual log total return: {actual_log_total:.4f}')
print(f'Sum of log returns:      {sum_of_log:.4f}')
print(f'Difference:              {abs(actual_log_total - sum_of_log):.6f}')
print()
print('They match! Log returns can be summed across time.')

# Convert between return types
print('=== Return Conversion ===')
print()

# Take first return as example
simple_r = simple_returns['SPY'].iloc[0]
log_r = log_returns['SPY'].iloc[0]

print(f'Original simple return:  {simple_r:.6f}')
print(f'Original log return:     {log_r:.6f}')
print()

# Convert simple to log
simple_to_log = np.log(1 + simple_r)
print(f'Simple → Log: ln(1 + {simple_r:.6f}) = {simple_to_log:.6f}')

# Convert log to simple
log_to_simple = np.exp(log_r) - 1
print(f'Log → Simple: e^{log_r:.6f} - 1 = {log_to_simple:.6f}')

# Visual comparison
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Distribution comparison
ax1 = axes[0]
ax1.hist(simple_returns['SPY'], bins=50, alpha=0.6, label='Simple', density=True)
ax1.hist(log_returns['SPY'], bins=50, alpha=0.6, label='Log', density=True)
ax1.set_xlabel('Return')
ax1.set_ylabel('Density')
ax1.set_title('Simple vs Log Returns Distribution')
ax1.legend()

# Scatter plot showing relationship
ax2 = axes[1]
ax2.scatter(simple_returns['SPY'], log_returns['SPY'], alpha=0.3, s=10)
ax2.plot([-0.15, 0.15], [-0.15, 0.15], 'r--', label='y = x')
ax2.set_xlabel('Simple Return')
ax2.set_ylabel('Log Return')
ax2.set_title('Simple vs Log Returns (nearly identical for small values)')
ax2.legend()

plt.tight_layout()
plt.show()

print('For small returns, simple ≈ log. For large returns, they diverge.')

# Compare different annualization methods
daily_mean = simple_returns['SPY'].mean()

print('=== Annualizing Daily Returns ===')
print()
print(f'Daily mean return: {daily_mean:.6f}')
print()

# Wrong way: simple multiplication
wrong_annual = daily_mean * 252

# Right way: compounding
right_annual = (1 + daily_mean) ** 252 - 1

print(f'WRONG (multiply by 252):  {wrong_annual:.2%}')
print(f'RIGHT (compound):         {right_annual:.2%}')
print(f'Difference:               {right_annual - wrong_annual:.2%}')
print()
print('For small returns, the difference is small. For larger returns, it matters!')

# Annualize volatility correctly
print('=== Annualized Volatility ===')
print()

for ticker in tickers:
    daily_vol = simple_returns[ticker].std()
    annual_vol = daily_vol * np.sqrt(252)
    print(f'{ticker}: Daily {daily_vol:.4f} → Annual {annual_vol:.2%}')

# Wrong vs right volatility annualization
daily_vol = simple_returns['SPY'].std()

wrong_vol = daily_vol * 252
right_vol = daily_vol * np.sqrt(252)

print('=== Volatility Annualization ===')
print()
print(f'Daily volatility:          {daily_vol:.4f}')
print()
print(f'WRONG (multiply by 252):   {wrong_vol:.2%}  ← Nonsensical!')
print(f'RIGHT (multiply by √252):  {right_vol:.2%}  ← Makes sense')
print()
print('Remember: √252 ≈ 15.87')

# Calculate Sharpe Ratio
risk_free_rate = 0.02  # Assume 2% annual risk-free rate

def calculate_sharpe(returns: pd.Series, risk_free_rate: float = 0.02) -> float:
    """Calculate annualized Sharpe ratio."""
    excess_returns = returns - (risk_free_rate / 252)
    return (excess_returns.mean() / excess_returns.std()) * np.sqrt(252)

print('=== Sharpe Ratios (Risk-Free Rate = 2%) ===')
print()

sharpe_ratios = {}
for ticker in tickers:
    sharpe = calculate_sharpe(simple_returns[ticker])
    sharpe_ratios[ticker] = sharpe
    quality = 'Excellent' if sharpe >= 1 else ('Good' if sharpe >= 0.5 else 'Poor')
    print(f'{ticker}: {sharpe:.3f} ({quality})')

def calculate_sortino(returns: pd.Series, risk_free_rate: float = 0.02, target: float = 0) -> float:
    """Calculate annualized Sortino ratio."""
    excess_returns = returns - (risk_free_rate / 252)
    
    # Downside deviation: only returns below target
    downside_returns = returns[returns < target]
    downside_std = np.sqrt(np.mean(downside_returns**2))
    
    return (excess_returns.mean() / downside_std) * np.sqrt(252)

print('=== Sortino Ratios (vs Sharpe) ===')
print()
print(f'{"Asset":<6} {"Sharpe":>10} {"Sortino":>10}')
print('-' * 28)

for ticker in tickers:
    sharpe = calculate_sharpe(simple_returns[ticker])
    sortino = calculate_sortino(simple_returns[ticker])
    print(f'{ticker:<6} {sharpe:>10.3f} {sortino:>10.3f}')

def calculate_max_drawdown(prices: pd.Series) -> float:
    """Calculate maximum drawdown from price series."""
    running_max = prices.cummax()
    drawdown = (prices - running_max) / running_max
    return drawdown.min()

def calculate_calmar(returns: pd.Series, prices: pd.Series) -> float:
    """Calculate Calmar ratio."""
    annual_return = (1 + returns.mean()) ** 252 - 1
    max_dd = calculate_max_drawdown(prices)
    return annual_return / abs(max_dd)

print('=== Maximum Drawdowns & Calmar Ratios ===')
print()
print(f'{"Asset":<6} {"Max Drawdown":>12} {"Annual Return":>14} {"Calmar":>10}')
print('-' * 46)

for ticker in tickers:
    max_dd = calculate_max_drawdown(prices[ticker])
    annual_ret = (1 + simple_returns[ticker].mean()) ** 252 - 1
    calmar = calculate_calmar(simple_returns[ticker], prices[ticker])
    print(f'{ticker:<6} {max_dd:>12.2%} {annual_ret:>14.2%} {calmar:>10.3f}')

# Visualize drawdowns
fig, axes = plt.subplots(2, 1, figsize=(14, 8))

# Normalize prices to start at 100
normalized = prices / prices.iloc[0] * 100

# Price chart
ax1 = axes[0]
for ticker in tickers:
    ax1.plot(normalized[ticker], label=ticker, linewidth=1.5)
ax1.set_ylabel('Normalized Price (100 = Start)')
ax1.set_title('Asset Prices Over Time')
ax1.legend(loc='upper left')

# Drawdown chart
ax2 = axes[1]
for ticker in tickers:
    running_max = prices[ticker].cummax()
    drawdown = (prices[ticker] - running_max) / running_max
    ax2.fill_between(drawdown.index, 0, drawdown, alpha=0.3, label=ticker)
ax2.set_ylabel('Drawdown')
ax2.set_title('Underwater Chart (Drawdowns)')
ax2.legend(loc='lower left')
ax2.yaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f'{x:.0%}'))

plt.tight_layout()
plt.show()

# Calculate Beta using SPY as the market
market_returns = simple_returns['SPY']

def calculate_beta(asset_returns: pd.Series, market_returns: pd.Series) -> float:
    """Calculate beta relative to market."""
    covariance = asset_returns.cov(market_returns)
    market_variance = market_returns.var()
    return covariance / market_variance

print('=== Beta Values (vs SPY) ===')
print()

for ticker in tickers:
    if ticker == 'SPY':
        beta = 1.0
    else:
        beta = calculate_beta(simple_returns[ticker], market_returns)
    
    interpretation = 'More risky' if beta > 1 else ('Less risky' if beta < 1 else 'Same as market')
    print(f'{ticker}: {beta:.3f} ({interpretation})')

def calculate_alpha(asset_returns: pd.Series, market_returns: pd.Series, 
                    risk_free_rate: float = 0.02) -> float:
    """Calculate annualized alpha using CAPM."""
    beta = calculate_beta(asset_returns, market_returns)
    
    # Annualized returns
    asset_annual = (1 + asset_returns.mean()) ** 252 - 1
    market_annual = (1 + market_returns.mean()) ** 252 - 1
    
    # Expected return based on CAPM
    expected = risk_free_rate + beta * (market_annual - risk_free_rate)
    
    # Alpha is the difference
    return asset_annual - expected

print('=== Alpha Values (Annualized) ===')
print()

for ticker in tickers:
    beta = calculate_beta(simple_returns[ticker], market_returns)
    alpha = calculate_alpha(simple_returns[ticker], market_returns)
    
    performance = 'Outperformed' if alpha > 0 else 'Underperformed'
    print(f'{ticker}: α = {alpha:>7.2%}, β = {beta:.3f} ({performance})')

def calculate_information_ratio(asset_returns: pd.Series, 
                                 benchmark_returns: pd.Series) -> float:
    """Calculate Information Ratio."""
    active_returns = asset_returns - benchmark_returns
    active_return_annual = active_returns.mean() * 252
    tracking_error = active_returns.std() * np.sqrt(252)
    return active_return_annual / tracking_error

print('=== Benchmark-Relative Metrics (vs SPY) ===')
print()
print(f'{"Asset":<6} {"Active Return":>14} {"Tracking Error":>14} {"Info Ratio":>12}')
print('-' * 50)

for ticker in tickers:
    if ticker == 'SPY':
        continue
        
    active_ret = (simple_returns[ticker].mean() - market_returns.mean()) * 252
    te = (simple_returns[ticker] - market_returns).std() * np.sqrt(252)
    ir = calculate_information_ratio(simple_returns[ticker], market_returns)
    
    print(f'{ticker:<6} {active_ret:>14.2%} {te:>14.2%} {ir:>12.3f}')